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O We identify the template organizing the chaotic dynamics 
^/"p f a vibrating string and show how to estimate topological pa- 
rameter values directly from an experimental time series. We 
x all these data processing techniques "topological time series 
analysis." The chaotic motions of a vibrating string are shown 
to be governed by a dissipative horseshoe map whose periodic 
CNorbit spectrum is, at the topological level, well described by 
' . ?L unimodal map of the interval. 

^■^ In addition, close topological agreement is demonstrated 
. between the experimental data and a synchronized model. 
I Thp empirical model of the experimental data allows us to cre- 
^^ate "synthetic data." Such synthetic data may enable a finer 
■"■^ — characterization of the system than that determined from the 
^seed data alone. 
>^ 
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O I- INTRODUCTION 

Our studies of the nonlinear and chaotic vibrations 
of a string are motivated by two intertwining concerns. 
First, we want to know strings. Or at least what we 
can know of strings via experiments showing chaotic mo- 
tions in a vibrating wire. Second, given an experimental 
times series from a low dimensional chaotic process, we 
would like to provide a quantitative characterization of 
the chaotic attractor — or at least those properties of the 
attractor which are possibly shared by a large class of 
similar chaotic systems. Such a characterization should 
be robust under the effects of noise and smooth parame- 
ter changes. Such a characterization is necessarily topo- 
logical. 

Interest in a string's motion has a long and illustrious 
history. In 1883 Kirchoff derived an inherently nonlinear 
differential equation which properly took into account the 
simultaneous longitudinal and transverse displacements 
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of a physical string [1]. In 1968 Narashima rederived a 
scaled version of Kirchoff 's equation [2]. Miles in 1984, 
beginning with Narashima's equation, derived and ana- 
lyzed a four dimensional model of the averaged oscilla- 
tions of a string and discovered a Hopf bifurcation [3]. In 
1989, after careful numerical studies, Johnson and Ba- 
jaj detected the presence of chaotic orbits in this same 
model [4]. Guided by these theoretical insights, exper- 
imental detection of these chaotic oscillations soon fol- 
lowed [5,6]. More complete historical references are found 
in several recent Ph. D. studies on the nonlinear vibra- 
tions of a string [7-9] . 

Of particular interest to our experimental studies are 
the theoretical examinations of low dimensional nonlin- 
ear string models by Miles [3], Bajaj and Johnson [10], 
Tufillaro [11], and O'Reilly [6,12]; and the experimental 
reports of chaotic vibrations in a string by Molteno and 
Tufillaro [5], Molteno [9], and O'Reilly and Holmes [6]. 
As remarked by most of these authors, chaotic motions 
in a string are not easy to observe. They exist only in 
a small parameter range and occur not at the forcing 
frequency, but rather in the slow oscillations of the am- 
plitude envelope. These amplitude modulations are usu- 
ally only a fraction of the overall transverse displacement. 
These facts may help to account for the rather late dis- 
covery of chaotic vibrations in such a common system. 

As recently advocated by several authors [13,14], the 
quantitative topological characterization of a low dimen- 
sional chaotic invariant set should proceed in two steps. 
First, the strange set needs to be assigned a good sym- 
bolic encoding. We address this problem by identify- 
ing the template organizing the stretching, twisting, and 
folding motion of the attractor in phase space [16]. Sec- 
ond, rules need to be developed specifying the presence or 
absence of subsets of symbol sequences when parameters 
are changed. The second problem is sometimes called 
pruning [15]. We present evidence that the "pruning 
problem" can be solved with predictions calculated by 



unimodal map theory — at least in the parameter regime 
considered and to within the resolution of the experi- 
mental measurements. Two types of evidence are used 
in support of these findings. The "symbol plane" is re- 
constructed for an experimental chaotic trajectory and it 
reveals that the system is well approximated by a vertical 
pruning front with no steps. Additionally, periodic orbits 
are extracted by the method of close recurrence [17,18] 
and their spectrum is ordered by unimodal theory. Since 
the data appears to be well described by unimodal the- 
ory, we also estimate the value of the kneading sequence, 
the "topological parameter" which fixes the spectrum of 
periodic orbits. 

In addition, we construct an empirical global model of 
the dynamics that produced the time series by determin- 
ing a vector which best fits the data [19]. The method 
used to fit the vector field is based on the minimum de- 
scription length (MDL) criterion of Rissanen [20]. Two 
distinct methods are used in an attempt to address the 
question of how closely the empirical model fits the data. 
First, we show that we can synchronize the model to the 
data. Synchronizing these models to the data set pro- 
vides the first piece of evidence that the empirical model 
is "close" to the fiow observed experimentally [27]. Sec- 
ond, we show that the experimental data and a time se- 
ries from the model share the same template and are close 
in terms of topological parameters. 

This paper is organized as follows. Section II describes 
the string experiment and examines data sets taken from 
a parameter regime which leads to the onset of a crisis. 
Section III briefiy describes how to construct a global 
empirical model from a data set. The Lyapunov spectra 
of the model and the data set are also calculated and it 
is suggested that calculating such system characteristics 
indirectly from "synthetic data" is a desirable and sen- 
sible procedure. Section IV describes how periodic orbit 
spectra are extracted by the method of close recurrence 
for data sets at two distinct parameter values. It also 
describes how braid invariants are used to identify the 
template organizing the periodic and chaotic orbits of the 
strange attractor. A topological analysis of data gener- 
ated from both the experiment and the empirical model 
is then presented. We call these data processing tech- 
niques "topological time series analysis," since they seek 
to ascertain topological form and topological parameter 
values directly from an experimental or simulated time 
series. Some concluding remarks are offered in Section V. 

Lastly, we would like to remark that our experimental 
results show good qualitative agreement with the numer- 
ical simulations on nonlinear and chaotic vibrations in 
strings by Bajaj and Johnson [10]. This correspondence 
is discussed more fully in Ref . [9] . 



II. EXPERIMENTAL SETUP AND DATA SET 
DESCRIPTION 

A schematic of the experimental apparatus used to ex- 
cite and detect forced vibrations in a thin wire is shown 
in Fig. 1. 
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FIG. 1. Schematic of experimental apparatus used to excite 
and detect forced vibrations in a thin wire. 

This apparatus is a considerably improved version of that 
described in Ref. [5]. A non-magnetic thoriated tung- 
sten wire, 0.15 mm in diameter, is fastened in a rigid 
mount. The wire sits in a permanent magnetic field cre- 
ated with rare earth ceramic magnets. An alternating 
current is passed through the wire to excite vibrations. 
Optoelectronic motion detectors are used to measure the 
transverse wire displacement simultaneously in two or- 
thogonal directions. The detectors consist of an infrared 
photodiode coupled to a matched phototransistor posi- 
tioned to detect the partial occlusion of the beam by the 
wire. 

Both the forcing current and detected signal are han- 
dled with a custom built controller/detector system de- 
signed around an Analog Devices ADSP 2105 Digital Sig- 
nal Processor (DSP). It is a 16-bit device with a pipelined 
Harvard Architecture capable of performing 30 million 
operations per second. Direct Digital Synthesis (DDS) 
is used to generate the signal for forcing the wire. This 
technique is used because fine control over the forcing 
frequency is required to characterize the constant ampli- 
tude response of the wire. This digital control and de- 
tection scheme also allows the development of real-time 
nonlinear signal processing and identification tools such 
as variable (in phase and amplitude) Poincare sections. 



Fast Fourier Transforms, and embedding plots with vari- 
able delays. All these real-time nonlinear signal process- 
ing tools greatly aid in the detection, identification, and 
optimization of the measurements reported here. 

Typical parameter values for the wire and appara- 
tus are: wire length, 0.07 m; mass per unit length, 
3.39 X 10-* kg m"*; density 2.1 x lO^* kg m"^; Young's 
modulus, 197778 x 10® N m~ ; magnetic field strength, 
0.2 Tesla; and a free frequency of vibration in the neigh- 
borhood of 1350 Hz. The damping of the system is mea- 
sured from the decay rate of free vibrations and can be 
modified by the application of a silicone coating to the 
wire [6]. 

The transverse amplitude displacement of the wire is 
usually sampled once each forcing cycle at a fixed phase 
(fd ~ 1300 Hz). The response frequency of the ampli- 
tude modulation is about 10 Hz. Thus, we have about 
130 samples per cycle — the system is over sampled. Such 
over sampling is necessary, however, in order to calculate 
accurately some of the braid (extracted periodic orbit) 
invariants described in the Section IV. Since the inter- 
esting dynamics occurs so slowly, it is possible to save 
large data sets directly to disk; data sets of 10® measure- 
ments are easy to obtain. The amplitude displacement 
measurements are usually low-noise (less than 1% of the 
RMS amplitude) and make good use of the wide dynamic 
range available. However, two measurement problems 
are present. First the system is subject to a parametric 
drift (typically manifested as a 1% variation of the RMS 
amplitude over 200 seconds) which we believe is due to 
a temperature instability in the drive system. Second, 
there are some systematic measurement errors due to 
phase jitter in the digital detecting electronics. Over- 
all, though, the measurements are of high quality: large, 
low noise, data sets making good use of the wide dynamic 
range made available from the 16-bit digitizers. 

After being recorded the raw data is "dynamically" 
cleaned by the nonlinear noise reduction technique de- 
scribed in Ref. [24]. The scalar data is embedded us- 
ing eleven dimensional delay coordinates. For each data 
point a small neighborhood is formed in which the two 
dimensional submanifold spanned by the attractor is ap- 
proximated by a linear subspace. Projections onto this 
subspace yield an improved time series. The procedure is 
iterated six times and curvature corrections are applied. 

In addition, dynamical cleaning provides an estimate 
of the noise level. For the data we have examined, the 
amplitude of the correction, which can be taken as an 
estimate for the amplitude of the noise in the data, is 
55 units (0.44%). We also estimate the amount of noise 
in the data independently by the method given in [25]. 
The functional form of the increase of the correlation in- 
tegral with embedding dimension is known for Gaussian 
measurement noise, so that the width of the Gaussian 
distribution can be determined by a simple function fit. 
We found that indeed the errors are compatible with a 
Gaussian distribution of width 65 units (0.5%). After 
nonlinear noise reduction the errors are no longer ex- 



pected to be Gaussian but the remaining noise level is 
found to be at most 25 units (0.2%). 

As the first step toward extracting a template from the 
time series, a three dimensional embedding of a chaotic 
trajectory is created out of the scalar amplitude measure- 
ments made by a single detector. This is accomplished 
via a time delayed embedding of the original scalar data 
set [21]. The offset for the delay is determined by the mu- 
tual information criterion [22]. Values for the offset range 
from 25 to 39 for the data sets examined, or, not unex- 
pectedly, about one quarter of the samples per chaotic 
cycle. A false nearest neighbor test also indicates that 
the time series is embeddable in three dimensions [21]. 
Figs. 2 shows three trajectories that are realized after a 
torus-doubling route to chaos [5]. Figs. 2(a) and 2(b) 
show chaotic attractors, and Fig. 2(c) shows a trajectory 
after the attractor has been extinguished by a crisis with 
a remote fixed point. 





FIG. 2. Three dimensional embedded time series from the 
string experiment at three different parameter values: (a) 
chaotic data set (a); chaotic data set (b); (c) crisis with a 
remote fixed point. The insets show the low dimensional char- 
acter of the next return maps constructed from a surface of 
section of the embedded attractors. 

These experimental results are in good qualitative agree- 
ment with the numerical simulations of Bajaj and John- 
son [10]. Fractal dimension calculations of these chaotic 
time series indicate a dimension of 2.1 [23]. Perhaps the 
strongest evidence, though, of the low dimensionality of 
these data sets is obtained by constructing a first return 
map of a two dimensional Poincare surface of section: 
these are shown in the insets of Figs. 2(a) and 2(b). 



III. A SYNCHRONIZED EMPIRICAL MODEL 

Our goal in this section is to create an empirical 
(global) model of the time series data. We do this by 
generating a vector field for a three dimensional system 
of ordinary differential equations whose dynamics "mim- 
ics" the dynamics of the experimental time series. The 
ability to "learn the evolution rules" only from a prion 
information is a powerful technique which opens the door 
to further applications such as short-term prediction and 
control. Also, as suggested here, the empirical model can 
be used to generate large, low noise, "synthetic" data 
sets. Such a synthetic data set, or the analytic knowl- 
edge supplied by the model, often permits a more refined 
calculation of a system's characteristics, such as its Lya- 
punov spectrum. However, before putting too much trust 
in the model, we must address the question of how faith- 
fully the empirical model represents the dynamics that 
generated the data. We approach this question by two 
distinct methods. In this section we show that the model 
can be synchronized to the data [26]. In the next section 
we compare topological properties extracted from the ex- 
perimental time series to the properties extracted from a 
time series generated by the model. 

The method used to generate an empirical model di- 
rectly from a data set is described in detail in a series 
of papers by Brown and co-workers [19]. The method 
assumes that the dynamics that produced the data vec- 
tors, y, can be written as a set of ordinary differential 
equations 



dy 
dt 



F(y). 



The vector field, F, is written as a expansion in terms of 
polynomials that are orthonormal on the attractor rep- 
resented by the experimental data 



F(y) = ^pm7rm(y). 



1=0 



These polynomials, tt*^'^ form a basis set in which to write 
a model. The expansion coefficients, p*^'^ are trained by 
modeling the system's dynamics as an implicit Adams 
integration scheme 



M 



y{t + r) = y{t) + r ^ aj^^F [y{t - {j - l)r)] 



j=0 



AM), 



where the a- ''s are the implicit Adams coefficients of 
order M. 

In the traditional numerical integration scheme the co- 
efficients, p'-'-''s, are viewed as being fixed. In contrast, 
in the modeling context, we view these coefficients as be- 
ing variable, their values being chosen to ensure a good 
correspondence between the experimental data and the 
empirical model. The model is trained by minimizing 



an MDL-type functional [20] which includes both a least 
squares minimization term and terms associated with the 
size and order of the vector field used to fit the data. The 
advantage of using an MDL-type functional is that an op- 
timal model, from the class of polynomial models, can be 
found. Heuristically, the parameterization chosen is op- 
timal in that it best fits the data with fewest number 
of terms. Higher order models may provide more accu- 
racy, but the cost of the additional terms in the MDL 
functional outweighs the profit gained by the increase in 
accuracy. 

To illustrate this modeling technique we consider the 
data set shown in Fig. 2(b). It consists of 128 x 10^ 
scalar 16-bit measurements of the vertical transverse 
string displacement. The forcing frequency in this par- 
ticular example is 1.384 kHz and the characteristic re- 
sponse frequency of the amplitude modulation is about 
9 Hz. The average mutual information and false near 
neighbor methods indicates that the optimal time delay 
and embedding dimension are T = 39 and d = 3 [21]. 
The model is trained by using a subsection of 10'* points 
of the entire data set. More details about the modeling 
of this particular data set can be found in Ref. [27]. 

As suggested in Ref. [19] after the model is constructed 
we then use the experimental time series to drive the em- 
pirical model. We find that it the experimental time se- 
ries is used to drive the model via dissipative coupling 
then the model synchronizes to the experimental time 
series. A detailed study of the quality of synchroniza- 
tion between the model and the data in the presence of 
additive noise and drift in the dynamics of the driving 
signal is presented in Ref. [27]. Synchronization between 
model and data provides the first piece of evidence that 
the empirical model and the experimental system are, in 
some respects, close. 

We have also calculated the spectrum of Lyapunov ex- 
ponents for the model as well as the data. Since it is the 
easier of the two calculations we first calculate the Lya- 
punov spectrum of the empirical model using a modified 
QR method [28-31]. We only report the two nonzero 
exponents of the three dimensional model. The results 
are shown in Figs. 3(a) and 3(b) where K represents the 
total number of Jacobians used for the calculations (the 
total evolution time is Kt where r = 0.05 is the normal- 
ized time between the data points). The figures indicate 
that for K > 5000 the values of Ai and A3 are essen- 
tially constant. We have used a small ensemble of ten 
different initial conditions taken from disparate regions 
of the attractor for our calculations. The error bars in 
the figures indicate the maximum and minimum values of 
the Lyapunov exponents calculated from this ensemble. 
The ensemble is small, however the error bars indicate 
that the variance of the calculated values of the Lya- 
punov exponents over the initial conditions is also small. 
Furthermore, if the initial conditions for our calculations 
are within the basin of attraction of the attractor then 
the Lyapunov exponents are independent of the initial 
conditions in the Kt -^ co limit [28,30]. 
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FIG. 3. Lyapunov spectrum for empirical model. 

The relative independence of the calculated values of 
the Lyapunov exponents for K > 5000 and different ini- 
tial conditions provide circumstantial evidence that our 
calculations have reached the asymptotic regime. We re- 
port the following values for the Lyapunov exponents of 
the empirical model: Ai = 4.31 x 10~^, and A3 = —0.576. 
These values are obtained by averaging over the ten ini- 
tial conditions for K = 15, 000. 

Calculating the full spectrum of Lyapunov exponents 
from an experimental data set is a notoriously difficult 
procedure. This is particularly true with regard to the 
negative exponents. In general, the best that one can 
hope for is to calculate a spectrum that is consistent with 
itself and whatever additional facts are known about the 
dynamics. The method we use to calculate the Lyapunov 
exponents from the experimental data is a minor varia- 
tion of the one reported in Ref. [32]. This method uses 
polynomials to map local neighborhoods on the attractor 
into their time evolved images. 

The first variation we employ involves using the mod- 
ified QR method reported in Refs. [30,31] instead of the 
one reported in Ref. [32]. The second modification is 
more complex and involves the data used to form the at- 
tractor. This second modification is required because the 



original data set is highly over sampled and needs to be 
"thinned out" before a reliable estimate of the spectrum 
can be obtained. The initial calculations use every other 
data vector to form the attractor. Next, every third data 
vector is used to form the attractor. Finally, every fourth 
data vector is used to form the attractor. 

The second modification also results in a change in 
the evolution time (we call it the step) associated with 
the polynomial maps from one local neighborhood to its 
image. We find that as the step increases the infiuence of 
experimental noise is reduced and consistent values are 
obtained for the Lyapunov exponents. In addition, the 
number of possible initial conditions increases. When 
we use every second data vector we have two different 
initial conditions, y(l) and y(2). Similarly when we use 
every third data vector we have three different initial 
conditions, etc. We use the different initial conditions as 
another consistency check for our calculations. For each 
case the number of Jacobians used in the calculation is 
K ~ N. 

In all cases the total number of vectors used to form 
the attractor is # ~ 10,000. Therefore, the size of the 
local neighborhoods used for the polynomial fitting is es- 
sentially the same for all tests. This insures consistency 
over our various tests. The value of # ~ 10, 000 is chosen 
on the ad hoc belief that larger values of N will result in 
local neighborhoods that are so small that their dynamics 
will be unduly infiuenced by noise in the data [33,34]. 

The results of these calculations are presented in 
Figs. 4(a) and 4(b). When the order of the fit is greater 
than two the calculated values of the Lyapunov expo- 
nents become independent of the order of the fit. The 
only exception occurs for A3 and the attractor obtained 
by using every other data vector, step = 2. Since these 
values are associated with the negative Lyapunov expo- 
nent, and are well outside the range of the other results, 
we will ignore them. As one would expect, there is some 
spread in the calculated values of the Lyapunov expo- 
nents depending on the initial condition used and the 
different evolution times (step = 2, 3, or 4). In general 
however, the spread is small for both the positive and 
the negative Lyapunov exponent. As a final matter we 
note that the total evolution time for our calculation is 
2Kt when we use every other data point, 3A'r when we 
use every third data point, etc. The independence of our 
results over step = 2, 3, and 4 indicates that K ~ 10, 000 
corresponds to the time asymptotic regime. The con- 
sistency of our results over changes in initial conditions, 
order of the fitting, and the evolution time leads us to 
conclude that the Lyapunov exponents calculated from 
the data set using this method are close to the true ex- 
ponents. 
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or "horseshoe" like map [15] is organizing the global dy- 
namics. Hence, one does not need the template to ascer- 
tain this topological information. However, the template 
contains an additional piece of topological information, 
namely, the global torsion, which can not be found by 
simply examining the first return map. 

To ascertain the form of the template we proceed in a 
straightforward way. First, we examine the data in the 
three dimensional embedding by rotating it, and viewing 
it from several viewpoints with the aid of a graphics pro- 
gram we have written [37]. This graphical visualization 
suggests that the sheeted structure shown in Fig. 5(a) 
properly portrays the stretching, twisting, and folding of 
phase space which is organizing the orbit structure. 




FIG. 4. Lyapunov spectrum from experimental data. 



To obtain numerical values of the Lyapunov exponents 
we first averaged the calculated values of A over the differ- 
ent initial conditions and then averaged that value over 
the order of the fit for all fits greater than 2. We find that 
Ai = 4.97 X 10-2 and A3 = -0.702. These values difl'er 
from those reported for the model by approximately 15% 
for the positive exponent and 20% for the negative expo- 
nent. 



IV. TOPOLOGICAL TIME SERIES ANALYSIS 
A. Template identification 

As the first step in a "topological time series analysis," 
we attempt to identify a "template" [16,35,36] organizing 
the periodic orbit structure of the fiow. A template is 
nothing more or less than a cartoon of the stretching 
and folding structure of the fiow written in a canonical 
form [38]. 

Before constructing the template, we note that the 
single-hump form of the first return map shown in the 
insets of Figs. 2(a) and (b) suggests that a "once-fold" 
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FIG. 5. Transformations on a horseshoe template: (a) 
Schematic of sheeted structure in the experimental data; 
(b) The template obtained from (a); (c) Horseshoe template 
equivalent to (a). 



The (negative) half-twist at the top of Fig. 5 is fol- 
lowed by a small (positive) fold shown at the bottom 
(compare with Fig. 5(b)). This small fold is responsi- 
ble for the stretching (causing sensitive dependence) and 
folding (causing the periodic and aperiodic recurrence) 
that produces the chaos. Both of these features of the 
flow are themselves organized by basic topological prop- 
erties of the flow, namely, the flow's fixed points and their 
homoclinic and heteroclinic connections. 

To arrive at the canonical form of the template we sep- 
arate the sheet shown in Fig. 5(a) into two branches by 
splitting along the trajectory of the fold point (in the 
language of Ref. [15], we are creating a symbolic parti- 
tion by considering the outermost "primary tangency"). 
Next we pull the small fold all the way to the bottom 
(thus going from a pruned to an unpruned system). The 
(negative) half twist in Fig. 5(a) results in a half twist 
in each branch, with both branches crossing as shown in 
the top of Fig. 5(b). The bottom of Fig. 5(b) shows how 
the small fold results in a horseshoe template in standard 
form [38]. Now we push all the branch twists in Fig. 5(b) 
to the top of the diagram and note that the twists in 
the left branch cancel, so that we arrive at the template 
shown in Fig. 5(c). We note that this is itself a horseshoe 
template with a global torsion of zero. Our construction 
thus shows that a standard horseshoe template with a 
global negative half twist is again a horseshoe template. 

In Subsection C we will extract periodic orbits from 
chaotic time series and discuss how each orbit is given a 
symbolic itinerary. Given these periodic orbits, we then 
attempt to verify (or falsify) the conjecture that a horse- 
shoe template with zero global torsion is organizing the 
dynamics. Specifically, we check two types of invariants: 
linking numbers of periodic orbits (a knot invariant), and 
the exponent sum (an invariant on positive braids) calcu- 
lated from the "natural" braid associated to each periodic 
orbit [38-40]. As noted by Hall [41], the exponent sum is 
a complete braid invariant for all horseshoe braids up to 
period eight — thus, this single invariant tells us the sym- 
bolic name (up to saddle-node pairs) of the low period 
orbits in a horseshoe without the need for an empirical 
symbolic partition. This is an effective procedure on or- 
bits all of whose crossings are clearly resolved. (These 
are all the orbits examined up to period eleven in the 
model, but only orbits up to period six in the experi- 
mental data.) We have found no discrepancies between 
the invariants predicted by a horseshoe template, those 
calculated from the data sets shown in Figs. 2, and data 
from the model constructed in Section IIL Some of these 
periodic orbits are shown in Fig. 6. These results pro- 
vide good evidence that — to within the resolution of the 
measurements and in the parameter regime considered — 
the orbits of the string and the model are governed by a 
horseshoe template with zero global torsion. 





FIG. 6. Some periodic orbits extracted from a chaotic time 
series from the string experiment. 

It is encouraging that the model agrees in this first 
topological check. Evidently, the analytic properties of 
the model class, along with the initial conditions provided 
by the experimental data, are sufficient to determine (at 
least approximately) the fixed points which are in the 
vicinity of the sampled fiow. 



B. Symbol Plane 

To get a topological "road map" of the data we use an 
empirical procedure to construct the symbol plane gen- 
erated by a single chaotic trajectory [15]. We construct 
two symbol planes, one for data from the experiment 
(Fig. 2(b)), consisting of 725 points in the return map, 
and one for the model, consisting of 4500 points in the re- 
turn map. To construct the symbol plane we first convert 
the chaotic trajectory to a symbolic string of O's and I's, 
depending on whether the orbit passes to the left or right 
of the maximum point of the first return map shown in 
the inset of Fig. 2(b). As shown in the next Subsection, 
this empirical symbolic partition is fine enough to distin- 
guish orbits at least up to period eleven. Now, since we 
know the topological form (the template) of the chaotic 
set, we can use kneading theory to determine the "well 



ordered" [15,42] symbol sequences needed to construct 
the symbolic plane. 

After conversion the scalar data set has become a sym- 
bol string of the form 

S = . . .S-3S-2S-IS0.S1S2S3 . . . 

where symbols to the left and right of sq are the past 
and future symbols respectively. The coordinates of the 
symbol plane are calculated from the well ordered past 
(cj) and future (hi) symbols as follows: 



<^) = J2^' ''i = J2' 



Dd2 




FIG. 7. Symbol planes generated by chaotic time series 
from a string: (a) Experimental data (b); Model (synthetic) 
data. 



D-1 



yi^) = Y.Vi 



8 = 



i + 1 

j=0 



mod 2. 



If s is an infinite symbol string generated by a chaotic 
orbit, then D is infinity in the above sums. However, 
since we are dealing with finite data sets, we approximate 
the symbol plane coordinates of a point by taking D = 
16. In this way we can use a finite symbol string from 
a chaotic trajectory to generate a sequence of points on 
the symbol plane. 

The resulting plots for both the experimental and 
model data are shown in Fig. 7. These plots suggest sev- 
eral nontrivial observations about the topological prop- 
erties of the fiows producing these data sets. First, 
it appears that a good approximation to the "pruning 
front" [15] is a single vertical line (a schematic for this 
pruning front is illustrated in Fig. 7), at least to within 
the resolution of our experimental measurements and 
data processing procedures. This indicates that, at the 
topological level, the symbolic dynamics of the system are 
well-described by unimodal map theory, with the prun- 
ing determined by a single topological parameter, the 
kneading invariant [15,42], which we can attempt to ap- 
proximate [43]. Second, the close similarity between the 
plot generated by experimental data, and model data, 
again provides striking evidence that the empirical model 
correctly captures the topological properties of the fiow. 
Moreover, we can estimate the closeness of this topolog- 
ical fit by comparing estimates of the kneading invari- 
ant from the experimental data and model data. Third, 
the simple form of these plots, and the conjectured (ap- 
proximate) pruning front, provides yet another powerful 
self-consistency check that the template (and hence the 
formula for well-ordering the symbol sequences) is iden- 
tified correctly. 



In the next section we extract periodic orbits from the 
chaotic time series and show how they can be used in 
this example to systematically approximate the vertical 
pruning front suggested by Fig. 7. 



C. Periodic orbit spectra 

To extract the (approximate) periodic orbits by the 
method of close recurrence we first convert the next im- 
pact map from the coordinate values directly into a sym- 
bol sequence. In this particular instance, we found that 
an adequate symbolic description (at least up to period 
eleven orbits, or approximately one part in 2^^) is ob- 
tained by choosing the maximum value of the next im- 
pact shown in the insets of in Fig. 2. Orbits passing to 
the left of the maximum are labeled '0', and those to 
the right are labeled '1'. Next we search this symbolic 
encoding for each and every periodic symbol string. Ev- 
ery time a periodic symbol string is found we calculate 
its (normalized) recurrence and then save the instance 
of the orbit with the best recurrence. The advantage of 
this procedure of periodic orbit extraction is that it is 
exhaustive. We search for every possible orbit up to a 
given period. In these studies we searched for all orbits 
of period one through eleven. 

The resulting spectrum of periodic orbits for both ex- 
perimental and model data sets is shown in Table 1. The 
orbits which are present in (the full shift) complete hy- 
perbolic system, and not present in the Table 1, are said 
to be pruned. Two topological invariants (the linking 
number and the exponent sum of positive braids) of the 
extracted orbits are calculated and compared with those 
of a horseshoe (see Table 1). There are no discrepan- 
cies on the orbits in which the crossings can be unam- 
biguously resolved. This indicates that, at least to this 
level of resolution, the template is correctly identified and 
the symbolic partition is adequate. Moreover, the tem- 
plate for the string experiment remains unchanged for 
the two different experimental parameter settings exam- 
ined, though the pruning is much more severe in data 
set (a) than it is in data set (c). The symbolic label 



(up to braid type) can also be determined — and used as 
yet another self-consistency check — by considering a sim- 
ple and easily computable braid invariants. As pointed 
out by Hall [41], the exponent sum (simply the sum of 
braid crossings in our example) is a complete invariant 
for horseshoe braids up to period eight. Also, an inspec- 
tion of the exponent sums as a function of period reveals 
that the exponent sum manages to distinguish most of 
the pseudo-Anosov horseshoe braids of periods nine, ten, 
and eleven as well [43]. The symbolics determined by the 
braid type, and that determined by the empirical parti- 
tion are consistent for all the orbits listed in Table 1. 
Our goal now is to predict as best as possible the pruned 
spectrum from the chaotic time series. 

TABLE I. Spectrum of low period orbits extracted from 
chaotic time series (all orbits with e < 0.03) are shown). 
Extracted orbits, their exponent sum, one-dimensional topo- 
logical entropy, and their (best) normalized recurrence are 
recorded. Data set (a) (experimental with 500 points in re- 
turn map); Data set (b) is the dynamically cleaned version 
of data set (c); Data set (c) (experimental with 725 points in 
return map); Data set (d) is from the empirical model of data 
set (c) — the "synthetic data" (4500 points in return map). 
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Before we discuss pruning, though, it is useful to con- 
sider the number of distinct periodic orbits extracted as a 
function of the number of points in the return map. This 
is shown for the model data in Table 2. As expected, the 
number of orbits that can be extracted increases with 
the number of points in the return map. More impor- 
tantly. Table 2 strongly suggests that using the method 
of close recurrence it is possible to obtain all the low pe- 
riod orbits embedded within the strange attractor. For 
instance, after 10® points are examined, we see that no 
new periodic orbits are found below period seven. Simi- 
larly, after examining 5 x 10® points we believe we have 
found all orbits up to period nine. These results also 
caution us when working with smaller samples such as 
those in data set (a) with 500 points in the return map, 
and data set (b) with 725 points in the return map — the 
extracted orbit spectrum is expected to miss orbits either 
because the orbit is pruned (it is not in the strange set) or 
because the sample of the strange set we are examining 
fails to provide a close enough coverage over the whole 
attractor. 

As suggested by the symbol plane diagram (Fig. 7), 
unimodal map theory should be sufficient to explain these 
periodic orbit spectra. To test this hypothesis, we first 
examine the periodic orbit spectra from the model data 
set and locate the maximal (rightmost) periodic point in 
terms of unimodal theory. It is a period five orbit with 
symbolic name 10110, and indicated by the label (up to 
saddle-node pairs) of s\. In unimodal theory, the s\ orbit 
forces [44] the orbits sfi, 4, SiD s\, sfi, si, s\-^, sfo, si, 
s\q, s\, s\q, s\, s\, s\, and s\. The notation identifies 
(up to saddle-node pairs) the period of each orbit in the 
subscript, and the value in terms of unimodal ordering 
(within each period) in the superscript. This theoreti- 
cally determined spectrum agrees with that found in Ta- 
ble 1. In fact, both partners of all the saddle-node pairs 
are present except for the period one orbit pair. In this 
case, though, it is known that the period one orbit with 
symbolic label '0' is not present in the strange attractor, 
though it is present in the fiow. 

TABLE IL Number of distinct periodic orbits of a given 
period extracted as a function of the sample size (x 10 ). 
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Moreover, we can now use the maximal periodic orbit to 
estimate the itinerary of the kneading sequence in the 
model data set, Km ~ 10110, or in two dimensions by the 
(vertical) pruning front specified by Oj^. 10110. The hori- 
zontal coordinate of the (maximal) point of this periodic 
orbit on the symbol plane is x(s) Ri 0.8485. An estimate 
of the kneading invariant also provides an estimate of the 
(one-dimensional) topological entropy [45,46]: hi(sl) k, 
0.4140. 

A similar analysis of the experimental data set (c) 
shows that the maximal periodic orbit in the unimodal 
ordering is s\i which forces Sg, s\i, s\, s\i, s\, s\i, s\q, 
s|, sfg, Sg, sjg, Sg, s\, s\, and s\. Up to period six, 
all the orbits forced by s\i are present. Nine orbits are 
missing between periods seven through eleven. But we 
believe that these higher period orbits are missing due to 
the small sample size. We use the maximal orbit to ap- 
proximate the itinerary of the kneading sequence in the 
experimental data set as kj Ri 10111101111, or in two 
dimensions by the (vertical) pruning front specified by 

0]^. 10111101111. The horizontal coordinate of the max- 
imal point of this periodic orbit on the symbol plane is 
x{s) Ri 0.8385. An estimate of the one-dimensional topo- 
logical entropy derived from the experimental data is: 
hi{s\^) f« 0.4032. 

Time series from the model and experiment are close 
in both their topological form and topological parameter 
value(s). The model, though, does appear to determine 
parameter values which over estimate the topological en- 
tropy. This difference in entropy could be due to several 
sources. First, uncertainties due to noise and parametric 
drift in the experimental data set place inherent limits on 
the modeling procedure's accuracy. Second, low-period 
orbits, perhaps of higher entropy than those extracted, 
may be present, but are not found because of the lim- 
ited sample size of the experimental data. Third, the 
estimation in the entropy can not be any finer than that 
allowed by the order of the periodic orbit approximation. 
As an examination of Table I indicates, the difference 
in entropy between the two periodic orbits considered 
(up to period eleven), which are close in the unimodal 
ordering, is roughly 0.01. This last point suggests that 
the estimated difference in entropy between the experi- 
mental data and model data is an upper bound on the 
difference. The entropy between the two process may, in 
fact, be closer than this estimate indicates. 



V. CONCLUDING REMARKS 

In this paper we create an empirical global model of 
a data set from a string experiment and show that the 
model captures the dynamics implicit within the data set 
in at least two significant respects: the model and data 
can be made to synchronize [19,26], and they share quan- 
tifiable common topological properties. We also suggest 
that it is a sensible and desirable procedure to use syn- 



thetic data from the model to characterize the system in 
ways which may not be readily accessible from the data 
set alone. In particular we discuss how to characterize 
the system at the topological level by its periodic orbit 
spectrum, and at the metric level by its Lyapunov spec- 
trum. More generally, we illustrate how to characterize 
a chaotic invariant set both by its general form (a tem- 
plate) and specific topological parameter values (knead- 
ing sequences) estimated from the spectrum of periodic 
orbits of the chaotic attractor, or from an empirically 
constructed "pruning front" [15]. 

Global, empirically constructed, analytic models open 
the door to a host of practical applications in systems 
analysis, identification, and control [19]. In particular, we 
emphasize how synthetic data sets permit the calculation 
of quantities requiring data with a larger sample size, 
or lower noise, than may be directly available from the 
experimental data alone [34]. In addition, we would also 
like to emphasize how a synchronized model of a system 
can be coupled directly to a data stream permitting a 
kind of model-based real-time dynamic filter [47]. 
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